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Abstract 

The fracture resistance of structures is optimised using the level-set method. Frac- 
ture resistance is assumed to be related to the elastic energy released by a crack 
propagating in a normal direction from parts of the boundary which are in tension, 
and is calculated using the virtual crack extension technique. The shape deriva- 
tive of the fracture-resistance objective function is derived. Two illustrative two- 
dimensional case studies are presented: a hole in a plate subjected to biaxial strain; 
and a bridge fixed at both ends subjected to a single load in which the compliance 
and fracture resistance are jointly optimised. The structures obtained have rounded 
corners and more material at places where they are in tension. Based on the results, 
we propose that fracture resistance may be modelled more easily but less directly 
by including a term proportional to surface area in the objective function, in con- 
junction with non-linear elasticity where the Young's modulus in tension is lower 
than in compression. 

Key words: Topology optimisation; Fracture; Level-set. 65K10 optimisation and 
variational techniques; 74R10 Brittle fracture; 74P10 optimisation of other 
properties; 74P15 Topological methods. 



1 Introduction 



Topology optimisation involves finding the geometry of an object which min- 
imises an objective function, with no constraints on the object's topology [1]. 
Broadly speaking, there are currently two main methods for studying such 
problems: the homogenisation method and its variants such as the well-established 



* Corresponding author. Telephone: +61 7 3365 3266. Fax: +61 7 3365 1477 

Email addresses: vchallis@maths.uq.edu.au (V J Challis), 
apr@maths.uq.edu.au (A P Roberts), awilkins@maths.uq.edu.au (A H 
Wilkins). 

Preprint submitted to Elsevier Science 6 February 2008 



SIMP method [2,3]; and, the level-set method [4]. In the former, the optimi- 
sation proceeds via a sequence of 'fictitious' materials which have properties 
unrealizable in nature until converging to the minimum containing only real 
materials [5-8]. In the latter, the process uses only real materials and moves 
boundaries to decrease the objective function at every iteration [9-12]. In 
its most naive implementation it typically converges slowly towards a local 
minimum, although recently some authors have proposed more sophisticated 
strategies to address these problems [13-24]. These strategies are not used 
here. 

Topology optimisation has been used to study many problems, and in partic- 
ular the elastic problem, where the objective function is the object's compli- 
ance or total elastic energy for a particular set of external loads and boundary 
constraints, has been well studied. Bends0e et al. [25] have prepared a recent 
review which discusses problems involving pressure loads, electromagnetic and 
acoustic wave propagation, fluid flow, articulated mechanisms and buckling. 

Level-set methods were first applied to structural optimisation in [9] and 
have since been applied to problems in two and three dimensions including 
compliance minimisation with single and multiple load cases [11,12,14,15, 
17-23,26-30] , multi-material problems [27,28,30,31], eigenfrequency prob- 
lems [10,19,29], gripping and clamping devices [11,15,18,19,26,27,31], and 
shape and image reconstruction [13,16]. 

The problem of optimising a material's fracture toughness is addressed in this 
paper. Our definition of fracture toughness involves consideration of cracks 
propagating in a normal direction from any point on the object's boundary 
into its interior. Thus, the level-set method is used since the objects of the 
intermediate iterations in the homogenisation/SIMP method have no well- 
defined "boundaries", being solid masses of fictitious materials. The use of 
this method also removes the problem of defining the fracture toughness of 
fictitious materials. 

The examples presented are two dimensional (2D). They are biphasic, consist- 
ing of "solid" which has a Young's modulus of unity and a Poisson's ratio of 
0.3, and "void" which has a Young's modulus of zero. 

Our approach closely follows Allaire et al. and for brevity we do not reproduce 
their methodology, rather we refer readers to the article [11]. 



2 



2 Brittle fracture and the objective function 

According to the classic work of Griffith [32] , a brittle material fractures if the 
elastic energy released upon crack propagation exceeds the material's "frac- 
ture surface energy" , which is a material-specific quantity parameterising the 
strength of the intermolecular bonds [33]. Given this quantity, the maximum 
allowed load for a particular situation may be obtained by calculating the 
elastic energy released for all possible crack propagation paths. A crack can 
nucleate and propagate from within a material, or it may nucleate from its 
boundary and propagate inwards. The latter situation only is considered here. 

The energy release rate for crack propagation may be calculated by finding the 
material's elastic energy for two configurations differing only by a small change 
in crack length. In practice, where the two energies are calculated numerically, 
this is computationally costly, and for the linear elastic case methods exist 
which only require one elastic relaxation, such as (for 2D) the J-integral [34], 
or equivalently, the virtual crack extension method [35,36]. 

Special cases of these methods may be interpreted as a particular shape deriva- 
tives as follows. Denote the domain occupied by the elastic material by Q. It 
has boundary dQ = U with outward unit normal n. The material is 
subjected to external tractions g on T N (Neumann conditions), body forces / 
in its interior, and its boundary has fixed displacement, u D , on T D (Dirichlet 
conditions). The displacement field that solves this linear elastic problem is 
denoted by Ui, and denote the energy density by W: 

W = \WiUjAijkiWkUi . 

Here and below the indices i,j,k,l = (1, . . . d) in d dimensions, summation over 
repeated indices is implied, the index n is reserved for the "normal direction" 
(e.g. x n = XiUi for any vector x), and Vj is the derivative, so the strain is 
\(ViUj + V jUi). The elastic tensor is denoted by A ijkl . The compliance is 



Consider deforming the material (coordinatised by x) with the diffeomorphism 
x — > x + 6{x) . 

Classical results [11,37-39] give the compliance change as a function of 9 as 




(1) 
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C (8) = 2 J 9-n (V n (g ■u) + Hg-u-W) 

T N 

+2 J 6-n (y n ( Uj - uf)A n3kl V kUl - W) . 



(2) 



In this formula V n denotes the derivative in the normal direction and H = V-n 
is the mean curvature of the boundary. 

To make the connection with the J-integral and virtual crack extension, con- 
sider the two dimensional (2D) case. Parameterise length along the boundary 
by y. (The boundary may consist of disjoint pieces, but that is irrelevant 
because of the following assumption.) Assume that 9 takes the form 



9 . n = ek (y) = < 



e(5 + y) for - 5 < y < 
e(S -y) for < y < 5 
otherwise 



on dQ. 



(3) 



Here e is a small number parameterising the "depth" of the boundary per- 
turbation, and S parameterises the "length" of the perturbation with 8 ^> e. 
The perturbation is shown in Fig. 1. The subscript on k indicates that the 




y = -5 y = y = 5 

Fig. 1. The perturbation of the boundary given by Eq. (3). 



perturbation is rooted at y — 0, and suppose that y = is part of T N . This 
is exactly a particular "virtual crack extension": in a finite element scheme 
the point y = is one node of the finite element mesh, ko corresponds to 
the virtual movement of this single node while keeping all other nodes fixed. 
Specialising Eq. (2) to this case gives 

C'{9) = J dy ek (y)(W - V n (g -u)-Hg-u) . (4) 



When the section of boundary — 5 < y < 5 is straight, H = 0, and this 
expression may be compared with Rice's J-integral evaluated along a contour 
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close to the point y — 0. 



In standard fracture situations where the material contains a well-defined crack 
of interest, numerical errors in the elasticity solution mean it is more accurate 
to move all the nodes which lie in the vicinity of a crack tip — k has nonzero 
support at all of these nodes — in order to calculate the energy release rate 
(or compliance change) for crack propagation in a certain direction. Similarly, 
for well-defined straight cracks, errors can be minimised by evaluating the J- 
integral around a contour that is distant from the crack tip. In the current 
application, however, generically there will be no well-defined cracks, and so 
the movement of just one node is considered, corresponding to a small crack 
initiating and propagating from the boundary in the normal direction. 

Motivated by these arguments, the objective function for fracture resistance 
considered here is 



The short-hand notation, G, has been introduced. The set S+ is all the points 
on Tjv where the surface tension is positive, so as to exclude points of high 
compressive energies which may not actually lead to fracture. Further com- 
ments on this point are made in Sec. 8. An explicit restriction to has been 
made, since the material is fixed at so that part of the boundary cannot be 
moved, so its fracture resistance cannot be optimised. When the exponent, q, 
is large, minimising J q corresponds to focussing on those few boundary points 
where the fracture resistance is low, while smaller q (say q — 1) considers 
points in S + more equally. The latter approach has been found to lead to 
more rapid convergence of the numerical algorithm, since the 'special points' 
in the former tend to move around during optimisation. 



3 The shape derivative of the objective function 

The shape derivative of the objective function is needed to implement the 
level-set scheme as detailed in Allaire et al. [11], since the velocity in the 
Hamilton- Jacobi (HJ) equation is the negative of the integrand of the shape 
derivative. 

To find the shape derivative of Eq. (5), the formal method discussed in [11] 
and attributed to Cea [40] is used. The idea is to find a Lagrangian, L, which 
generates the equations of motion for the field(s) and reduces to the objective 
function when the fields satisfy their equations of motion ("on shell"). Because 




where 




(5) 
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of the former property, the derivative <5L/<5fields = on shell, so the shape 
derivative dL/dfl can be computed at fixed fields. 

The elasticity problem is defined through the equations of motion 



Vj(i4 i jfciV fc «i) = -/» in ^, 
A njM V k ui = gj on T N , 

Mi = uf on T D . (6) 

Written in terms of stress, a, the first is simply VjCr^- = — and the second 
is a n j = gj, so / is a body force and g a boundary traction. 

There are many Lagrangians which generate these equations of motion, and 
the following contains just one auxialiary field M which is a Lagrange multi- 
plier: 

L M = J WiM 3 A i]kl W k ui - J MJi - J M l3i 

-J [MjAjfcjVfcUi + (uj - uf )A njkl V k M l ] . (7) 

Varying Mj gives the correct equations of motion for u; and 
Lm = , 



on shell. 

Therefore, the Lagrangian 
L = J q + Lm , 



satisfies the requirements given above. Varying u in this Lagrangian gives the 
equations of motion for the auxiliary field: 



Vj(i4y fc ,V fc M,)=0 in n, 

A njkl V k M t = V* {k p qG q - l A ljkl V kUl ) + k p qG q - 1 (y n g 3 + H 9j ) on T 
Mi = on F D . 

The derivative V* acts in the tangent plane only: 

V* = Vi - UiVn . 
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In deriving these equations, the boundary conditions on u have been used, and 
it is assumed that k p has zero support on dT^, so that there is no boundary 
term when integrating by parts with the operator V*. These equations there- 
fore define an elasticity problem with zero displacement on r fl and tractions 
applied to the boundary around the points p G S+. 

Using the standard formulae found in [11], and evaluating on-shell (so that 
u = u D and M = on T D , for instance), the shape derivative of the objective 
function is 

J' q (0) = J (Vn (k p G q ) + Hk p G q ) 9-n + J V \M 3 A idk tf \u x 9 ■ n 
s + an 

-J (V n (g-M)+Hg-M)0-n 

- J (V n MjA njk iV k ui + VnUjAyfcjVfcM,) 9 ■ n. (9) 



4 Compliance minimisation and the volume or area constraint 



In addition to fracture resistance, compliance minimisation is also considered 
in the examples below. The objective function used is 

J = (l-X)C + XJ q , (10) 



where A is a fixed weighting factor. Both C and J q are minimised when the 
material has infinite volume (3D) or area (2D). To avoid this, a maximum 
volume (3D) or area (2D), Umax, is introduced, and the problem to be solved 
is 

Minimise J subject to J 1 < V max . (11) 

n 



It is common in the literature to add a term proportional to the volume (or 
area) to the objective function, so the problem becomes: minimise a linear 
combination of J and volume (or area). Unfortunately this does not allow di- 
rect comparison of structures with different A since altering A will also produce 
different volume fractions. The problem stated in Eq. (11) does allow a direct 
comparison. 
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5 Finite element implementation 



The finite element method is implemented using bi-linear (in 2D) square el- 
ements, as described in [41-43]. The original (u) elastic problem is solved, 
and G at each boundary node is calculated to find the boundary tractions 
for the auxiliary (M) problem. Four geometric arrangements are possible, as 
illustrated in Fig. 2. In the case of Fig. 2(a), the integral becomes 

5 

J k p G q = J dy k (G q 22 - G q 21 ) x=0 + jdy k (G q 12 - G q u ) x=0 , (12) 

-S 



where the subscripts 11, etc, correspond to the values in the elements defined 
in Fig. 2. In calculating these G values, H = since the boundary is straight, 
and in the case illustrated Gn = = G 2 \. In principle, if there is a boundary 
traction, acting at the central node then its action gets split over the 
elements 12 and 22 at y = (however, in the practical situations considered 
here and elsewhere the boundary never moves away from such a point so the 
point may be simply excluded from S+). 

In the case of Fig. 2(b) and (c) the normal is taken to be (— 1, — 1)/a/2 to 
ascertain whether the point is in S + . Algorithmically it is convenient to use 
the sum of two orthogonal virtual crack extensions as indicated in the figure, 
each with a similar expression to Eq. (12). In the case of Fig. 2(d), the nor- 
mal direction is not well-defined. In this case, the two normal virtual crack 
extensions are tried (each can be decomposed into two orthogonal parts as in 
(b) and (c)), and the maximum of the two results (which will be numerically 
similar) is chosen. 



(a) 



ii 



21 



\ 12 



22 



(b) 



(c) 



e 
x/2 




V2 



Fig. 2. Four possibilities for a node (black dot) on the boundary, (a): there are 
two neighbouring solid elements and two void elements, (b): three solid elements, 
(c): one solid element, (d): two diagonally-opposite solid elements. Virtual crack 
extensions for each are shown by dotted lines: (b), (c) and (d) can be decomposed 
into a linear combination of two orthogonal extensions. The surrounding elements 
are labelled 11, 12, 21 and 22 as shown in (a). 
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The elasticity problem is then solved for M. The surface tractions are applied 
at each point p G S + and its two neighbouring nodes on the boundary. A 
finite difference scheme is adopted to find the tangential divergence in Eq. (8): 
k v qG q ~ x Gij is evaluated at points half way between the nodes (y = ±5/2 with 
x = in Fig. 2(a)) and the appropriate differences taken. Cases Fig. 2(b), (c) 
and (d) follow similarly for instance, in case (b) 



x=-S/2,y=0 



- k qG q - 



-1 ~ a 2jj 



V2 



x=0,y=-S/2/ 



since the tangential stress is (aij — a 2 j)/V2 and the distance between the 
points (-5/2, 0) and (0, -5/2) is 5/^2. 

After the elasticity problem for M has been solved, the shape derivative of 
J q can be found by straightforward evaluation. The only nontrivial term in 
Eq (9) is V n (k p G q ). For Fig. 2(a), this is evaluated by by taking the finite 
difference of Eq. (12) at x — 1 (using a Young's modulus of zero for the 12 
and 22 elements in the virtual crack extension calculation) with that at x — 0. 
Similar constructions are used for the other 3 cases in Fig. 2. 

The void material has zero Young's modulus and thus the integrand of the 
shape derivative is zero for void elements. Consequently, the velocity, v, of 
the HJ evolution is zero too, and so the boundaries of the object will never 
move! Allaire et al. [11] promote the use of a so-called ersatz material, where 
the voids have a small nonzero Young's modulus, allowing the velocity to 
be extended into the void elements. This is not suitable here, since the first 
and third terms of Eq. (9) are only defined on the boundary. Therefore, a 
smoothing operation is performed as follows. The first and third terms are 
defined at the nodal points, so their contribution is spread equally to the four 
neighbouring elements. The second is evaluated at the midpoints of each of 
the solid elements, and then spread to the neighbours by the correlation: 

Vncw(i, j) = (2v(i,j) + v(i - I J) +v(i + I J) +v(i,j - 1) +v(i,j + l))/6 



This smoothing of the boundary values into the void allows the boundary to 
move during HJ evolution. 

As found by other authors, the objective function typically decreases with 
each HJ iteration, only increasing when the solid volume fraction exceeds V max . 
When that occurs, the material is eroded away using a constant velocity in the 
HJ evolution until the volume fraction is 0.99Vmax- The behaviour of J q is not 
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as controlled as in the compliance-minimisation case because sometimes the 
movement of one part of the boundary causes G to increase markedly some- 
where else. This may be seen in Fig. 6. Because our finite-element program 
is parallelised and hence rather quick, we evolve very slowly: each iteration 
involves an elastic relaxation and movement of the boundary by only one 
element. 



6 Case study 1: the hole 



The fracture resistance of a plate containing an arbitrarily shaped central hole 
must be maximised under the conditions that the hole area fraction is 1/8 and 
the plate is under uniform biaxial strain. This problem may be solved using 
A = 1 and V mSuX = 7/8. Fig. 3 shows one initial condition, and the geometry at 
the minimum, which has a circular hole. Because this is the expected result, 
this case study illustrates that the algorithm is indeed maximising fracture 
resistance. 




(b) 




Fig. 3. (a): An initial configuration with J2 = 973. (b) The minimum at J2 = 177. 
This example has 200 x 200 pixels. 



7 Case study 2: the bridge 



The 'bridge' problem is specified in Fig. 4 which also shows its solution for 
A = 0, VVax = 17.1%, and the length variable being twice the height. Similar 
problems were studied in [1,11,12,14,15,18,21,22]. The compliance of the 
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structure shown in Fig. 4(b) is C = 26.5. This topology and geometry is also 
obtained using the SIMP method [1]. The simulations were performed on a 
200 x 100 mesh. A large variety of initial configurations were used in order to 
address the problem of converging to a local minimum, rather than a global 
one. 

The minimum structure at A = has J\ = 0.5 — 50 times less than C - 
so to weight fracture resistance meaningfully a value of A ~ 1 must be used. 
The following study considers the extreme case of A = 1 — 10~ 5 . Less extreme 
values produce geometries that interpolate between the one shown here and 
the A = case. In order to keep the structure within the rectangular domain 
of interest, it is assumed that any parts of the bridge on the bottom boundary 
are infinitely fracture resistant, otherwise the optimisation results in material 
enveloping the point of application of the load. 

Figure 5 shows the optimal geometry, and Fig. 6 shows the objective function. 
This geometry has C = 34.3, 3 X = 0.222 and so J = 0.222. The following 
points may be noticed. 

(1) The parts of the structure which are in tension capture material from 
those in compression. 

(2) Corners are typically more rounded. 

(3) In the case of the bridge, the optimal topology is unchanged. In other 
cases that have been investigated (such as the cantilever), this is some- 
times not true since a thin strut in compression can disappear entirely. 

(4) G is fairly constant over the parts of the boundary that are in tension, as 
shown in Fig. 7. This suggests that very similar geometries are obtained 
for q > 1 (although the value of A must also be scaled to obtain exactly 
the same geometry). 

These points are supported by further simulations with different q and A: 
because the structures obtained are actually quite similar we resist reproducing 
further graphics for brevity. 



8 Discussion 

In this paper fracture resistance has been defined through a virtual crack ex- 
tension in a normal direction on parts of the boundary that are in tension. 
Because the normal stresses are typically close to zero on the boundary (even 
considering the few places where g ^ 0, discretisation effects and numerical 
errors which all yield nonzero normal stresses), the objective function for frac- 
ture resistance is close to an integral of the energy density over parts of the 
boundary, S+, which are in tension. 
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If this restriction to S+ was removed, we propose that the the objective func- 
tion for fracture resistance may be well approximated by a surface integral of 
the energy density, and therefore, when used in linear combination with the ob- 
jective function for compliance, which is an integral of the energy density over 
the volume of the material, may be further replaced by an objective function 
proportional to the surface area (in 3D) or perimeter (in 2D). This provides a 
physical motivation for using a surface area or perimeter constraint/objective 
function, which are typically employed to penalise micro-fine structures some- 
times encountered in structural optimisation problems. 

The problem is more complicated with the restriction to S+, as material is 
captured by the parts in tension. Those parts do not necessarily decrease 
significantly in surface area, but because they are more massive will have sig- 
nificantly less surface energy. Our experience is that in most cases the optimal 
topology does not change compared with the compliance-only situation, al- 
though the geometry may be significantly altered, as can be seen in Fig. 5. 
Similarly to above, we propose that this case may be modelled less directly, 
but probably more easily by including a term proportional to the surface area 
or perimeter in the objective function, in conjunction with using a non-linear 
elasticity law where the Young's modulus in tension is lower than in compres- 
sion. 
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• i-H 
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Fig. 4. (a): The bridge problem has two points with fixed zero displacement and a 
single load, (b): The solution for width = 2 x height, V max = 17.1% and A = 0. 




Fig. 5. Left: The optimal geometry for the bridge problem for A = 1 — 10 6 . Right 
a comparison of the left-hand structure and the optimal geometry for A = 0. 
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Fig. 6. The objective function at each iteration for the bridge case study. Each 
iteration involves either: adding two elements (symmetrically about the symmetry 
axis) to the solid structure; or, eroding the structure so that the solid volume fraction 
is 0.99V max . The latter causes large jumps while the former typically decreases J, as 
discussed in the text. The initial configuration had topology and geometry similar 
to that of the compliance-optimised solution of Fig. 4(b). 



Fig. 7. Greyscaling indicating the value of G for each boundary node in the optimal 
solution shown in Fig. 5. Dark regions, such as those close to the supporting points 
indicate high G, lighter regions indicate lower G. 
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